#include <iostream>
#include <sstream>
#include <cstring>
#include <vector>
#include <queue>
#include <unordered_map>
#include <algorithm>
#include <random>
#include <array>
#include "getopt.h"
#include <sys/stat.h>
#include "config-fnt.h"
#include "graph.priv.h"
#include "stream.h"
#include "public/utils.h"

namespace std {
	template<typename T, std::size_t N> struct hash<std::array<T, N>> {
		std::size_t operator()(const std::array<T, N>& k) const {
			std::size_t h{0};
			for(std::size_t i=0; i<N; i++) {
				h^=std::hash<T>{}(k[i]);
			}
			return h;
		}
	};
}

/* Options */
static const char* opts=":o:d:r:hv";
static const struct option opts_long[]=
{
	{"output", 1, nullptr, 'o'},
	{"dist", 1, nullptr, 'd'},
	{"ratio", 1, nullptr, 'r'},
	{"help", 0, nullptr, 'h'},
	{"version", 0, nullptr, 'v'},
	{nullptr, 0, nullptr, 0}
};

/* Print version information and exit */
void version(const char* name) {
	std::cout<<name<<" (" PACKAGE_NAME_LONG ") " PACKAGE_VERSION "\n";
	std::cout<<"Copyright (C) 2015-2016 " PACKAGE_ORG "\n";
	std::cout<<"License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>\n\n";
	std::cout<<"Written by: " PACKAGE_AUTHORS "\n";
	exit(0);
}
// //
/* Print usage information and exit */
void usage(const char* name, int ecode) {
	if(ecode) {
		std::cerr<<"Try '"<<name<<" --help' for more information.\n";
		exit(ecode);
	}
	std::cout<<"Usage: "<<name<<" [OPTION]... FNT_FILE0 FNT_FILE1\n";
	std::cout<<
		"Calculate FNT file differences.\n\n"
		"Options:\n"
		"  -o, --output=FILE [OPTIONAL] Save joined files to FILE. By default,\n"
		"                       save to standard output.\n"
		"  -r, --ratio=RAT   [OPTIONAL] Distance-to-length ratio for checking"
	  	"                       edge duplication.\n"
		"  -d, --dist=DIST   [OPTIONAL] Distance for checking edge duplication.\n\n"
		"  -h, --help        [OPTIONAL] Print this help, and do nothing else.\n"
		"  -v, --version     [OPTIONAL] Print version information and exit.\n\n"
		"Report bugs to <" PACKAGE_BUGREPORT ">\n"
		PACKAGE_NAME_LONG " home page: <" PACKAGE_URL ">\n";
	exit(ecode);
}


static int cmp_point(const Point& a, const Point& b) {
	if(a._x<b._x) return -1;
	if(a._x>b._x) return 1;
	if(a._y<b._y) return -1;
	if(a._y>b._y) return 1;
	if(a._z<b._z) return -1;
	if(a._z>b._z) return 1;

	if(a._r<b._r)
		throwError("not possible");
	if(a._r>b._r)
		throwError("not possible");
	if(a.m<b.m)
		throwError("not possible");
	if(a.m>b.m)
		throwError("not possible");
	return 0;
}

static int cmp_edge(Edge a, Edge b) {
	auto epa=EdgePriv::get(a);
	auto epb=EdgePriv::get(b);
	if(epa->type<epb->type)
		return -1;
	if(epa->type>epb->type)
		return 1;
	if(epa->points.size()<epb->points.size())
		return -1;
	if(epa->points.size()>epb->points.size())
		return 1;
	auto n=epa->points.size();
	for(size_t i=0; i<n; i++) {
		auto r=cmp_point(epa->points[i], epb->points[i]);
		if(r<0)
			return -1;
		if(r>0)
			return 1;
	}
	return 0;
}

static int cmp_edge_idx(Edge a, Edge b) {
	auto epa=EdgePriv::get(a);
	auto epb=EdgePriv::get(b);
	if(epa->index<epb->index)
		return -1;
	if(epa->index>epb->index)
		return 1;
	return 0;
}

typedef std::pair<Edge, bool> EdgeEnd;
static int cmp_edge_end(const EdgeEnd& a, const EdgeEnd& b) {
	auto r=cmp_edge_idx(a.first, b.first);
	if(r<0)
		return -1;
	if(r>0)
		return 1;
	if(!a.second && b.second)
		return -1;
	if(a.second && !b.second)
		return 1;
	return 0;
}

static int cmp_vert(Vertex a, Vertex b) {
	auto vpa=VertexPriv::get(a);
	auto vpb=VertexPriv::get(b);
	if(!vpa->finished && vpb->finished)
		return -1;
	if(vpa->finished && !vpb->finished)
		return 1;

	if(vpa->neighbors.size()<vpb->neighbors.size())
		return -1;
	if(vpa->neighbors.size()>vpb->neighbors.size())
		return 1;
	auto n=vpa->neighbors.size();
	for(size_t i=0; i<n; i++) {
		auto r=cmp_edge_end(vpa->neighbors[i], vpb->neighbors[i]);
		if(r<0)
			return -1;
		if(r>0)
			return 1;
	}
	return 0;
}

static int cmp_nrn(Tree a, Tree b) {
	auto vpa=VertexPriv::get(a.root());
	auto vpb=VertexPriv::get(b.root());
	if(vpa->index<vpb->index)
		return -1;
	if(vpa->index>vpb->index)
		return 1;
	return 0;
}

struct Component {
	std::vector<size_t> edges;
	std::vector<size_t> verts;
	std::vector<size_t> nrns;
};

void reorder_component(Component* comp, GraphPriv* grp) {
	std::sort(comp->edges.begin(), comp->edges.end(), [grp](size_t a, size_t b)->bool { return cmp_edge(grp->edges[a], grp->edges[b])<0; });
	for(size_t i=0; i<comp->edges.size(); i++) {
		auto ep=EdgePriv::get(grp->edges[comp->edges[i]]);
		ep->index=i;
	}
	for(size_t i=0; i<comp->verts.size(); i++) {
		auto vp=VertexPriv::get(grp->vertices[comp->verts[i]]);
		std::sort(vp->neighbors.begin(), vp->neighbors.end(), [](const EdgeEnd& a, const EdgeEnd& b)->bool { return cmp_edge_end(a, b)<0; });
	}
	std::sort(comp->verts.begin(), comp->verts.end(), [grp](size_t a, size_t b)->bool { return cmp_vert(grp->vertices[a], grp->vertices[b])<0; });
	for(size_t i=0; i<comp->verts.size(); i++) {
		auto ep=VertexPriv::get(grp->vertices[comp->verts[i]]);
		ep->index=i;
	}
	std::sort(comp->nrns.begin(), comp->nrns.end(), [grp](size_t a, size_t b)->bool { return cmp_nrn(grp->trees[a], grp->trees[b])<0; });
}

bool eq_component(Component* a, Component* b, GraphPriv* grp) {
	if(a->edges.size()!=b->edges.size())
		return false;
	for(size_t i=0; i<a->edges.size(); i++) {
		auto r=cmp_edge(grp->edges[a->edges[i]], grp->edges[b->edges[i]]);
		if(r!=0)
			return false;
	}
	if(a->verts.size()!=b->verts.size())
		return false;
	for(size_t i=0; i<a->verts.size(); i++) {
		auto r=cmp_vert(grp->vertices[a->verts[i]], grp->vertices[b->verts[i]]);
		if(r!=0)
			return false;
	}
	if(a->nrns.size()!=b->nrns.size())
		return false;
	for(size_t i=0; i<a->nrns.size(); i++) {
		auto r=cmp_nrn(grp->trees[a->nrns[i]], grp->trees[b->nrns[i]]);
		if(r!=0)
			return false;
	}
	return true;
}


std::vector<Component> find_components(Graph gr) {
	std::vector<Component> comps;
	auto n=gr.vertices().size();
	std::vector<bool> edg_dirty(gr.edges().size(), false);
	std::vector<bool> vert_dirty(n, false);
	std::vector<bool> nrn_dirty(gr.trees().size(), false);
	while(true) {
		size_t seed=SIZE_MAX;
		for(size_t i=0; i<n; i++) {
			if(!vert_dirty[i]) {
				seed=i;
				break;
			}
		}
		if(seed==SIZE_MAX)
			break;
		Component comp{{}, {}, {}};
		std::queue<size_t> que{};
		que.push(seed);
		while(!que.empty()) {
			auto i=que.front();
			que.pop();
			auto v=gr.vertices()[i];
			if(!vert_dirty[i]) {
				comp.verts.emplace_back(i);
				for(auto& p: v.neighbors()) {
					auto v1=p.second?p.first.leftVertex():p.first.rightVertex();
					que.push(v1.index());
				}
				auto nrn=v.tree();
				if(nrn && nrn.root()==v) {
				}
				vert_dirty[i]=true;
			}
		}
		for(auto i: comp.verts) {
			auto v=gr.vertices()[i];
			for(auto& p: v.neighbors()) {
				auto j=p.first.index();
				if(!edg_dirty[j]) {
					comp.edges.push_back(j);
					edg_dirty[j]=true;
				}
			}
		}
		for(auto i: comp.verts) {
			auto nrn=gr.vertices()[i].tree();
			if(nrn) {
				auto j=nrn.index();
				if(!nrn_dirty[j]) {
					comp.nrns.push_back(j);
					nrn_dirty[j]=true;
				}
			}
		}
		comps.push_back(std::move(comp));
	}
	for(auto v: edg_dirty) {
		if(!v)
			throwError("Untouched edge");
	}
	for(auto v: vert_dirty) {
		if(!v)
			throwError("Untouched vertex");
	}
	for(auto v: nrn_dirty) {
		if(!v)
			throwError("Untouched neuron");
	}
	return comps;
}

inline double calcDist(const Point& p0, const Point& p1, const Point& p2) {
	double p0x=p0.x();
	double p0y=p0.y();
	double p0z=p0.z();
	double p1x=p1.x();
	double p1y=p1.y();
	double p1z=p1.z();
	double p2x=p2.x();
	double p2y=p2.y();
	double p2z=p2.z();
	double ux=p1x-p2x;
	double uy=p1y-p2y;
	double uz=p1z-p2z;
	double vx=p2x-p0x;
	double vy=p2y-p0y;
	double vz=p2z-p0z;
	double a=ux*ux+uy*uy+uz*uz;
	double b=ux*vx+uy*vy+uz*vz;
	double c=vx*vx+vy*vy+vz*vz;
	double t=-b/a;
	if(t<0) t=0;
	if(t>1) t=1;
	return sqrt(a*t*t+2*b*t+c);
}


template<bool D>
bool is_dup(Edge e0, Edge e1, double r, double r2, bool alien) {
	auto n0=e0.points().size();
	auto n1=e1.points().size();
	auto idx0=[n0](size_t k) ->size_t { return D?n0-1-k:k; };
	auto& pts0=e0.points();
	auto& pts1=e1.points();
	if(alien) {
		if(cmp_point(pts0[idx0(0)], pts1[0])!=0
				&& cmp_point(pts0[idx0(n0-1)], pts1[n1-1])!=0)
			return false;
	} else {
		if(cmp_point(pts0[idx0(0)], pts1[0])!=0
				|| cmp_point(pts0[idx0(n0-1)], pts1[n1-1])!=0)
			return false;
	}
	double maxdist=0.0;
	size_t k0, k1;
	for(k0=0, k1=0; k0<n0; k0++) {
		auto dist=calcDist(pts0[idx0(k0)], pts1[k1], pts1[k1+1]);
		for(auto k11=k1+1; k11+1<n1 && k11<k1+8; k11++) {
			auto d=calcDist(pts0[idx0(k0)], pts1[k11], pts1[k11+1]);
			if(d<dist) {
				dist=d;
				k1=k11;
			}
		}
		if(dist>maxdist)
			maxdist=dist;
	}
	double len=0.0;
	for(k1=1; k1<n1; k1++)
		len+=pts1[k1-1].distTo(pts1[k1]);
	
	if(maxdist<=r && maxdist<=r2*len)
		return true;
	std::cerr<<"dist: "<<maxdist<<std::endl;
	std::cerr<<"rat: "<<maxdist/len<<std::endl;
	return false;
}
int main(int argc, char* argv[]) {
	const char* output=nullptr;
	double rad{0.0};
	double rat{0.0};

	// Parse options
	int opt;
	int prev_optind=optind;
	while((opt=getopt_long(argc, argv, opts, opts_long, nullptr))!=-1) {
		switch(opt) {
			case '?':
				std::cerr<<argv[0]<<": unrecognized option '";
				std::cerr<<argv[prev_optind]<<"'\n";
				usage(argv[0], -1);
			case ':':
				std::cerr<<argv[0]<<": option '"<<argv[prev_optind];
				std::cerr<<"' requires an argument\n";
				usage(argv[0], -1);
			case 'h':
				usage(argv[0], 0);
			case 'v':
				version(argv[0]);
			case 'o':
				output=optarg;
				break;
			case 'd':
				{
					char* eptr;
					rad=strtod(optarg, &eptr);
					if(*eptr!='\x00') {
						std::cerr<<"Failed to parse dist\n";
						usage(argv[0], -1);
					}
				}
				break;
			case 'r':
				{
					char* eptr;
					rat=strtod(optarg, &eptr);
					if(*eptr!='\x00') {
						std::cerr<<"Failed to parse ratio\n";
						usage(argv[0], -1);
					}
				}
				break;
			default:
				std::cerr<<"Unexpected error.\n";
				usage(argv[0], -1);
		}
		prev_optind=optind;
	}

	if(!output) {
		std::cerr<<"Where to save output?\n";
		usage(argv[0], -1);
	}
	if(optind+2!=argc) {
		std::cerr<<"Need 2 input FNT files.\n";
		usage(argv[0], -1);
	}

	std::vector<std::string> header;
	auto gr=GraphPriv::alloc();
	auto gp=GraphPriv::get(gr);

	struct MapVal {
		std::array<bool, 2> hits{false, false};
		bool isvert{false};
		double r{0.0};
		int16_t m{0};
		int dup{0};
	};
	using MapKey=std::array<int32_t, 3>;

	std::unordered_map<MapKey, MapVal> points;
	std::unordered_map<EdgePriv*, bool> edge2alien{};
	for(int j=0; j<2; j++) {
		auto gr0=GraphPriv::alloc();
		auto gp0=GraphPriv::get(gr0);
		auto hdr=gp0->loadFnt(argv[optind+j]);
		if(j==0)
			header=std::move(hdr);

		for(auto e: gr0.edges()) {
			auto& pts=e.points();
			for(size_t i=0; i<pts.size(); i++) {
				auto& p=pts[i];
				MapKey key{p._x, p._y, p._z};
				auto ins=points.emplace(key, MapVal{});

				if(i==0 || i==pts.size()-1) {
					ins.first->second.isvert=true;
				}

				ins.first->second.hits[j]=true;
				ins.first->second.r+=p.r();
				ins.first->second.m|=p.m;
				ins.first->second.dup+=1;
			}
			gp->edges.push_back(e);
			if(j)
				edge2alien[EdgePriv::get(e)]=true;
		}
		for(auto v: gr0.vertices()) {
			gp->vertices.push_back(v);
		}
		for(auto v: gr0.trees()) {
			gp->trees.push_back(v);
		}
		GraphPriv::free(gr0);
	}

	do {
		std::vector<Edge> edges{};
		do {
			bool split{false};
			for(auto e: gr.edges()) {
				auto ep=EdgePriv::get(e);
				std::vector<size_t> brks{};
				MapVal* prev_val{nullptr};
				auto& pts=ep->points;
				for(size_t i=0; i<pts.size(); i++) {
					auto& p=pts[i];
					auto& v=points.at(MapKey{p._x, p._y, p._z});
					p.r(v.r/v.dup);
					p.m=v.m;
					size_t brk;
					MapVal* brk_val{nullptr};
					if(i>0 && prev_val->hits!=v.hits) {
						if(prev_val->hits==std::array<bool, 2>{true, true}) {
							brk_val=prev_val;
							brk=i-1;
						} else if(v.hits==std::array<bool, 2>{true, true}) {
							brk_val=&v;
							brk=i;
						} else {
							throwError("not possible");
						}
					} else if(v.isvert && i!=0 && i+1!=pts.size()) {
						brk_val=&v;
						brk=i;
					}
					if(brk_val) {
						if(brk==0) {
							//auto vp=VertexPriv::get(e.leftVertex());
							// XXX
							//vp->finished=false;
						} else if(brk+1==pts.size()) {
							//auto vp=VertexPriv::get(e.rightVertex());
							// XXX
							//vp->finished=false;
						} else {
							if(brks.empty()) {
								brk_val->isvert=true;
								brks.push_back(brk);
							} else {
								if(brks.back()<brk) {
									brk_val->isvert=true;
									brks.push_back(brk);
								}
							}
						}
					}
					prev_val=&v;
				}
				if(brks.empty()) {
					edges.push_back(e);
					continue;
				}
				split=true;
				std::vector<Vertex> verts{};
				bool alien=edge2alien[EdgePriv::get(e)];
				for(size_t k=0; k<=brks.size(); k++) {
					Vertex v0, v1;
					size_t a, b;
					if(k==0) {
						v0=e.leftVertex();
						a=0;
					} else {
						v0=verts[k-1];
						a=brks[k-1];
					}
					if(k==brks.size()) {
						v1=e.rightVertex();
						b=pts.size()-1;
					} else {
						v1=VertexPriv::alloc();
						auto v1p=VertexPriv::get(v1);
						//v1p->finished=false;
						// XXX
						v1p->finished=true;
						verts.push_back(v1);
						b=brks[k];
					}
					auto e0=EdgePriv::alloc();
					auto e0p=EdgePriv::get(e0);
					e0p->type=e.type();
					for(size_t j=a; j<=b; j++)
						e0p->points.push_back(pts[j]);
					e0p->leftVertex=v0;
					e0p->rightVertex=v1;
					if(alien)
						edge2alien[e0p]=true;
					edges.push_back(e0);
				}
				for(auto v: verts)
					gp->vertices.push_back(v);
			}
			if(!split)
				break;
			std::swap(gp->edges, edges);
			edges.clear();
		} while(true);

#if 0
		std::swap(gp->edges, edges);
#else
		std::unordered_map<EdgePriv*, bool> edge2dup{};
		gp->edges.clear();
		for(size_t i=0; i<edges.size(); i++) {
			bool dup{false};
			auto e0=edges[i];
			bool alien=edge2alien[EdgePriv::get(e0)];
			for(size_t j=0; j<i; j++) {
				auto e1=edges[j];
				dup=(is_dup<true>(e0, e1, rad, rat, alien) && is_dup<true>(e1, e0, rad, rat, alien))
					|| (is_dup<false>(e0, e1, rad, rat, alien) && is_dup<false>(e1, e0, rad, rat, alien));
				if(dup) {
					edge2dup[EdgePriv::get(e1)]=true;
					bool leftSoma{false}, rightSoma{false};
					for(auto t: gr.trees()) {
						if(t.root()==e0.leftVertex())
							leftSoma=true;
						if(t.root()==e0.rightVertex())
							rightSoma=true;
					}
					bool preserved{false};
					if(leftSoma || rightSoma) {
						if(e0.points().size()<64) {
							edge2dup[EdgePriv::get(e0)]=true;
							gp->edges.push_back(e0);
							preserved=true;
						} else {
							if(leftSoma) {
								auto e2=EdgePriv::alloc();
								auto e2p=EdgePriv::get(e2);
								e2p->type=e0.type();
								for(size_t ii=0; ii<16; ii++)
									e2p->points.push_back(e0.points()[ii]);
								e2p->leftVertex=e0.leftVertex();
								auto vv=VertexPriv::alloc();
								auto vvp=VertexPriv::get(vv);
								vvp->finished=false;
								e2p->rightVertex=vv;
								gp->vertices.push_back(vv);
								edge2dup[EdgePriv::get(e0)]=true;
								gp->edges.push_back(e2);
							}
							if(rightSoma) {
								auto e2=EdgePriv::alloc();
								auto e2p=EdgePriv::get(e2);
								e2p->type=e0.type();
								for(size_t ii=e0.points().size()-16; ii<e0.points().size(); ii++)
									e2p->points.push_back(e0.points()[ii]);
								e2p->rightVertex=e0.rightVertex();
								auto vv=VertexPriv::alloc();
								auto vvp=VertexPriv::get(vv);
								vvp->finished=false;
								e2p->leftVertex=vv;
								gp->vertices.push_back(vv);
								edge2dup[EdgePriv::get(e0)]=true;
								gp->edges.push_back(e2);
							}
						}
					}
					if(!preserved) {
						auto e1p=EdgePriv::get(e1);
						if(e1p->type<e0.type())
							e1p->type=e0.type();
						auto v0p=VertexPriv::get(e0.leftVertex());
						v0p->finished=false;
						auto v1p=VertexPriv::get(e0.rightVertex());
						v1p->finished=false;
					}
					break;
				}
			}
			if(!dup) {
				edge2dup[EdgePriv::get(e0)]=false;
				gp->edges.push_back(e0);
			}
		}
#endif

		for(auto v: gp->vertices) {
			auto vp=VertexPriv::get(v);
			vp->neighbors.clear();
		}

		for(auto e: gp->edges) {
			auto v0=e.leftVertex();
			auto v0p=VertexPriv::get(v0);
			v0p->neighbors.emplace_back(e, false);
			auto v1=e.rightVertex();
			auto v1p=VertexPriv::get(v1);
			v1p->neighbors.emplace_back(e, true);
		}
		{
			std::vector<Vertex> verts{};
			for(auto v: gp->vertices) {
				if(!v.neighbors().empty())
					verts.push_back(v);
			}
			std::swap(gp->vertices, verts);
		}
		for(auto v: gp->vertices) {
			auto vp=VertexPriv::get(v);
			bool has0{false}, has1{false};
			for(auto& p: vp->neighbors) {
				if(edge2dup[EdgePriv::get(p.first)]) {
					has1=true;
				} else {
					has0=true;
				}
			}
			// XXX
			if(has1 && has0)
				vp->finished=false;
		}
	} while(false);

	{
#if 0
		std::default_random_engine rng{};
		std::shuffle(gp->edges.begin(), gp->edges.end(), rng);
		std::shuffle(gp->vertices.begin(), gp->vertices.end(), rng);
		std::shuffle(gp->trees.begin(), gp->trees.end(), rng);
#endif

		for(auto e: gp->edges) {
			auto ep=EdgePriv::get(e);
			ep->type=0;
		}
		gp->updateModel();
	}

	if(output) {
		gp->saveFnt(output, header);
	} else {
		gp->save(std::cout);
	}
	GraphPriv::free(gr);

	return 0;
}

